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SUMMARY 

Algorithms have been developed using Monte Carlo techniques to determine the thermal 
network parameters necessary to perform a finite difference analysis on Combinatorial Solid 
Geometry (CSG) models. Orbital and laser fluxes as well as internal heat generation are modeled to 
facilitate satellite modeling. The results of the thermal calculations are used to model die infra-red 
(IR) images of targets and assess target vulnerability. Sample analyses and validation are 
presented which demonstrate code products. 


INTRODUCTION 

CSG models of targets have been used for many years in performing various phenomenologic 
analyses such as nuclear particle transport, kinetic energy weapon effects, survivability, susceptibility, 
and lethality studies. Many CSG vehicle target description databases are under construction and 
include highly detailed, three dimensional solid geometry models of red and blue, high value strategic 
and tactical targets. Ideally, these data sets include all of the information necessary to perform 
detailed thermal analysis for thermal performance simulation. Applications of thermal analysis on 
these data sets include thermal design analysis, laser weapon effects survivability analysis, and IR 
image synthesis for sensor and seeker performance of target susceptibility, acquisition and tracking, 
simulation of system IR imaging, and automatic target recognition algorithm training. 

SINDA and other network methods for thermal analysis algorithms have been successfully 
employed to model many complex thermophysical systems. Their success may be partially credited 
to the simplicity of modeling very complicated thermal control, radiative transfer and nonlinear 
thermal problems. However, these codes are not usually linked directly to geometric solid geometry 
models, and therefore the burden of constructing the input data has been placed on the code user. 
Often, the only cross correlation for temperature predictions is a sketch or drawing and a node 
number penciled in by hand. 

At the Phillips Laboratory and Battelle Memorial Institute, algorithms have recently been 
developed to marry the two technologies, allowing high resolution thermal analysis of these CSG 
target sets for such studies. The algorithms include three dimensional geometry discretization, 
thermal network parameter computation including interelemental radiative transfer computation. 
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orbital and external thermal environment simulation, transient source computation including internal 
heat generations and thermal control. The computed data is then reformatted to a SINDA input file 
format for transient or steady state thermal predictions. The resultant thermal profiles are explicitly 
linked to the geometric models. This correlation increases the usefulness of the data by allowing IR 
image and signature predictions from the total set of geometry, thermophysical, thermoradiative and 
temperature data. Although originally created for imagery analysis, the algorithms were modified to 
include time and space resolved laser fluxes. The thermal results of the laser loads on a target can be 
used in vulnerability analysis. All of the algorithms are combined in a suite of codes collectively 
named AutoSINDA. 


THERMAL SOLVER CHOSEN 

SINDA was chosen as the thermal solver for AutoSINDA for a variety of reasons. Its 
application is flexible, able to model time and temperature dependent parameters, convection, 
conduction, internal and external radiation, heat pipes, thermal control systems, and phase changes. 
SINDA also allows programmable access through user defined variables and has many solution 
schemes which allow both transient and steady state solutions. In addition, SINDA has been used and 
tested for many years and has a heritage as a thermal analysis code in the satellite community. The 
SINDA code has a broad base of users within the community from which to draw expertise. 


ALGORITHM DEVELOPMENT 
Overview 

An overview of the AutoSINDA codes is presented in Figure 1. The model is first discretized. 
The discretized model has the same format as the original model and has lost no geometric fidelity. 
Next the thermal network parameters are calculated using a Monte Carlo technique. Orbital loadings 
are applied using only first incident absorption of solar and earth fluxes. Diffuse fluxes such as earth 
shine is modeled as a distributed impressed flux. The result is an element/load table as a function of 
time. AH thermal network parameters and loads are reformatted for SINDA input. SINDA is then 
used for transient thermal analysis. There is a one to one correspondence between the discretized 
e ements in the solid model and the SINDA nodes. The temperatures correlated to the discretized 
geometry are used to determine laser effects or the resulting IR image of the target. 


Assumptions 

"^ e AutoSINDA model is limited by the following assumptions: 1) all materials are opaque 
2) all materials are perfect diffuse or perfect specular, and 3) only first incident external radiative flux 
absorption is modeled. 


CSG Geometric Description 

A combinatorial solid geometry model is described by Boolean combinations of solid 
primitives. The lowest hierarchy unit of these models is called the "element" (not related to element 
as in Finite Element). Each element should be composed of a contiguous solid material and have a 
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single surface property or coating. An example of a single element would be a base plate for an 
antenna or a panel for a structural housing. 

Each element is a built by a combination of solid primitives. This entails complex intersections 

and unions of void and solid ellipsoids, elliptical cones and polyhedra. Note that in general these 

combinatorial intersections of volumes allows very complex solid descriptions in a compact format. 

The target model is comprised of many elements assembled in a series of hierarchial sets to 

describe the complete system model. Many of these models are quite detailed, modeled down to the 

electronic card level. 


Algorithm Basic Processes 

In order to perform thermal analyses on the CSG descriptions, five basic processes are 
executed on the geometry model. These are: 

Process 1: Model Discretization 

Process 2: Thermal Network Parameter Computation 

Process 3: Environmental Simulation 

Process 4: Application of Environmental Flux & Internal Heat Generations 

Process 5: Data Reformatting to SINDA Input Deck 

Model Discretization 

For each element in the CSG model, a bounding box is found using the minimum and 
maximum x, y, and z coordinates for the element. The number of divisions in each direction which 
will produce a length less than or equal to the maximum mesh size, for that element is determined 

by, 

n r “ + 1 (1) 

where r is the axis being divided and "int" takes the truncated integer of the expression. Each axis is 
divided into n, equal segments. A plane going through each division is cut perpendicular to the axis. 
Planes which cut through the x, y, and z axis intersect to form new boxes whose sides are at most 
long. Each box created that overlaps with any part of the original element is "ANDed" with the 
original element using boolean algebra. This creates a new element in the discretized model. This 
algorithm discretizes all elements in this "rectangular" manner and may create "strange" or very small 
new elements. If an element from the original model has a convex surface, this discretization scheme 
is undesirably susceptible to creating elements that are in two discontiguous pieces. In this case, 
increasing the discrete resolution will reduce this anomaly. 

Thermal Network Parameter Computation 

A Monte Carlo technique is employed to determine volumes, surface areas, conductive 
resistances, and radiative resistances. For each element, I, a sphere is found which completely 
surrounds the element. The sphere is found by determining the minimum and maximum coordinates 
for each primitive in the element (an element can be composed of several primitives which are 
"ANDed" and "ORed" together). The maximum coordinate for the element is determined by the 
maximum coordinate of all element primitives that are "Ored" and the minimum maximum coordinate 
for all the primitives that are "ANDed". The minimum coordinate for the element is found similarly. 
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Note that this technique is very simply and efficiently evaluated but does not necessarily find the 
smallest sphere containing the element, which may be extremely complex. 

The radius of this sphere is determined by, 

U _ ~ + fomc ~ * few ~ ^mfa) (2) 
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The center of the sphere is located at the average coordinate. 
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The cross-sectional area through the center of the sphere is given by, 

A — - 
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Rays are shot in an inward random direction from a random location on the surface of the 
surrounding sphere. A total of n^ are shot. The number of times a ray hits the surface of the 
element, n^,, is recorded. Every time a ray hits the surface of the element, there is one entrance and 
one exit point. For every hit, the distance traveled through the element from the entrance to the exit 
point, dL, is recorded. The total length traveled through the element is determined by summing the 
distances traveled through the element from each hit. 
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The volume of the element is determined by, 
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The surface area of the element is determined by, 
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For every entrance and exit point that a ray found, a ray is shot from that location on the 
surface of the element in the outward normal direction. The ray is tracked until it either hits another 
element or escapes to space. If the ray hits another element, k, the distance traveled to the element k 
is determined. If the distance traveled is less than a tolerance distance for conduction, the point 
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where the ray originated is located on a conductive surface. The number of times a point on a 
conductive area is found from the element I to element k, is collected in n^,^. If the origination 
point is located on a conductive surface, a ray is shot from that location in the inward normal 
direction until it exits the element. The distance traveled to the exit of element I is recorded in the 
array d k . When the ray exits element I, a check is made to see if it conducts to the same element, k, 
at that location. If it does, d k is divided by two. The distance to the conductive center from element 
I to k will be determined by half the average value of d*. The reason d k is divided by two if the ray 
conducts to element k on both sides of element I, is that instead of heat flowing through element I, 
heat will flow in the inward/outward direction from the element. 

For every entrance and exit point found that is not located on a conductive area, the energy 
that is. emitted from the element to all the other elements and space is found. A ray is shot in a 
random outward direction. The ray is given an initial energy weight, w, of: 
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The ray is tracked until it either hits another element or escapes to space. If the ray escapes to space, 
the remaining weight, «, is added to and the ray is terminated. If the ray hits an element, k, 
the energy deposited in E,.* is determined by, 

E lk " E i* + e *“ <9) 


The weight of the ray is reduced by the energy deposited in the element it hit, 

o ■ u - e k a 


( 10 ) 


The ray continues in a random direction from the point it hit on element k until an exit criteria is met. 
Exit criteria are in the form of a maximum number of reflections made by the ray and fraction of 
original weight left in the ray. 

The thermal network from element I to all other elements and space is determined. For every 
element k that element I conducts to, the conductive area between I and k is found by. 
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The average length, l^v.,!-*. from the "center" of element I that is normal to the conductive area, 
Acood.i-k is found by. 
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The non-linear radiative conductance (inverse of resistance) from element I to element k is contained 
in E,_*. 

Upon completion of this algorithm for every element, the arrays from i to j and j to i are 
combined. Since the conductive area between element i and j is computed separately for element i 
and element j, they must be combined. A weighted average between the two areas is taken. They 
are weighted by the cross-sectional area of the other’s surrounding cylinder. 
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The conductive resistance from element i to element j can then be determined by, 
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where is the thermal conductivity. The "resistive" conductor between element i and element j is 
found by adding the two resistors in series from the elements to the conductive area, 

+ Rjm (^ 


The conductive heat flow between element i and j is then found by taking the inverse of the 
resistance. 

To find the radiative heat flow between element i and j, a weighted average between and 
Ej^j is taken. They are weighted the same way as the areas, 

E H “ w i e h + ( 16 ) 


Environmental Simulation 

Environmental simulation is accomplished using a standard orbital code which reads 'the orbital 
elements and epoch of the satellite and writes the satellite ephemeris to a data file. This information 
includes the direct solar flux and direction, the earth shine flux, the solar albedo flux and the 
orientation of the satellite with respect to the earth. 


Application of External Fluxes & Internal Heat Generations 

The direct sun, laser flux, and earth albedo are modeled as collimated loads. Although the 
earth albedo is by no means collimated, this model is implemented for computational simplicity. This 
model does conserve the energy of the earth albedo, primarily making the gross approximation in the 
directional distribution. 

Earthshine is modeled in a distributed fashion. The distributed earthshine is modeled as 
originating from a flat earth disk, broken into 13 sections of equal area. Each area radiates one 
thirteenth of the total earth flux on the space object. The disk radius is defined as, 

r d - sinQ (17) 


where 
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0 - aumfa? + 2 a 


( 18 ) 


and 


altitude of satellite 


( 19 ) 


and r e is the IR apparent radius of the earth, 

r t - 6378.2 fan + 30 hn(IR apparent atmosphere) (20) 

The disk is broken into an inner disk, a four-piece middle annulus, and an eight-piece outer annulus. 
The area of each of the thirteen sections is, 


* section 


13 


( 21 ) 


The radius of the inner circle is. 



( 22 ) 


The outer radius of the middle annulus is, 



The middle annulus is sectioned at 0, 90, 180, and 270 degrees. The outer annulus is sectioned at 0, 
45, 90, 135, 180, 225, 270, and 315 degrees. 

The effective directions for each section to impinge the satellite is modeled deterministically as 
coming from a point at the center of each sector. The center is defined by the area average radius for 
the annulus sectors and the center of the circle for the inner circle sector. Each direction is assigned 
one thirteenth of the total earthshine flux. 

The flux and direction of all external loads are applied to the solid model by ray tracing. 
Currently only first-incidence actions of externally applied loads are modeled. The absorbed flux 
going into each element is numerically evaluated by ray tracing the entire model in a semi-systematic 
uniform grid. The grid is defined by the square which surrounds the circle defined by the model- 
enveloping sphere. The enveloping sphere is defined as the sphere centered on the models bounding 
rectangular box with a radius extending to the bounding box’s comer. Ray tracing is accomplished 
through a user-specified mesh size. 

The radiant heat transfer into each element from loads other than laser loads is calculated by 
multiplying the number of hits on the element, the area each ray represents, the spectral (or banded) 
absorptivity of each element, and the uniform incident flux of the load. Since laser loads do not have 
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uniform flux, the laser model employs rays with an energy equal to its representative area times the 
radial variation of flux according to the gaussian laser profile, 


9 



(24) 


where q,, is the peak flux (central intensity) of the laser, r is the radial distance from the ray being 
traced to the aim point of the target, and a is the Gaussian deviation parameter for characterizing the 
spot size. 

Internal Heat Generations 


Internal loadings can be specified as a constant heat generation or an on/off thermal control 
heat generation. The internal heats are assigned to the geometry model’s original elements before 
discretization and are then distributed to the discrete elements based on mass fraction. 

Data Reformat to SINDA Input Deck 

Upon completion of the four previous processes, all necessary data to perform thermal analysis 
now exists. The data is formatted to a SINDA compatible input deck and subsequently submitted to 
the SINDA process. This reformatting operation allows the use of both temperature dependent or 
independent material properties. The option to cyclically repeat a single orbit is also allowed. Any 
of the SINDA transient solution schemes may be used as well as the steady state solvers. 


SAMPLE CASES 
Conductive Model Validation 
Case 1: Solid Cube with a Zero Temperature Boundary Condition 

Two samples cases were created to test the validity of the conduction method, the first sample 
case is a box with sides length 5 cm. The temperature of the box is initially 100 K when a zero 
temperature is imposed at the boundary for times greater than 0. In order to apply the above 
algorithms, a thin shell .1 cm thick is modeled around the box. The analytical solution for a box with 
sides of 5 and 5.1 cm is compared to the solution predicted by the algor ithms f or differe nt mesh sizes 
in Table 1. In the table, n is the number of divisions in the x, y, and z directions. 

As the mesh size gets finer (n gets larger), the algorithms approach the analytical solution with 
the exception of the finest discretization. The final column in Table 1 shows a divergent behavior as 
the number of divisions was increased to 17 in each direction. This is probably caused by numerical 
round-off errors. 

Case 2: Sphere with a Zero Temperature Boundary Condition 


The second sample case is a sphere with radius 5 cm initially set to 100 K when a zero 
temperature is imposed at the boundary for times greater than 0. Again, in order to apply the above 
algorithms, a thin shell 0.1 cm thick is modeled around the sphere. The analytical solution for the 
sphere with radii of 5 and 5.1 cm is compared to the solution predicted by the algorithms for different 
mesh sizes in Table 2. 
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Table I Comparison of the temperature at the center of the box between the analytical solution 
and the numerical approximation 


time 

(s) 

anal 

(5cm) 

anal 

(5.1cm) 

CO 

tl 

c 

n=6 

n=9 

n= 1 1 

n= 17 









1 

69.8 

72.0 

88.2 

79.7 

74.2 

71.6 

73.8 

2 

26.9 

29.0 

33.2 

21.7 

29.8 

28.8 

30.4 

3 

9.8 

11.0 

0.0 

m 

11.2 

10.9 

11.7 

4 

3.5 

4.1 

0.0 

2.3 

4.2 

4.1 

4.5 

5 

1.3 

1.6 

0.0 

0.7 

1.6 

1.5 

1.7 

6 

0.5 

0.6 

0.0 

0.2 

0.6 

0.6 

0.7 

7 

0.2 

0.2 

0.0 

0.1 

0.2 

0.2 

0.3 

8 

0.1 

0.1 

0.0 

0.0 

0.1 

0.1 

0.1 


Radiative Transfer Validation 

In order to validate the radiative transfer model, a simple test geometry was constructed using 
the solid geometry modeling software. A 6 sided hollow box was created 3 x 4 x 5 cm in internal 
dimensions x, y and z respectively. Note that the radiative transfer mechanism accounts for both 
internal and external radiative heat transfer, so it was necessary to ensure that this box was 
constructed with no external view factors. Two methods were used to compute the radiative heat 
transfer coefficients. 

The first method involved using a 6 facet analytical model for 6 simple view factors. This 
method assumes uniform distribution of all irradiance and radiosity. The method of Gebhart was used 
to find the radiative interchange factors and the matrix. 

The second method involved discretizing each face of the box and numerically computing view 
factors for each small element within the faces. Upon determining the hundreds of view factors, the 
method of Gebhart was used to compute the radiative transfer factors. The large matrix of radiative 
interchange factors then had to be summed over each face using the appropriate summation principles. 
Note that the AutoSINDA algorithms inherently model the second case more accurately as it does not 
assume the uniform irradiance and radiosity distribution. 

The two distinct methods were chosen to allow complete comparison, knowing that the second 
method models the radiative transfer more accurately. Table 3 summarizes the model of the 6 sided 
box and its 300K broad band radiative properties. 

Case 1: 6 Node Analytical Radiative Transfer 

The internal radiative interchange factors for this geometry, K^, were computed using a 
simple 6 node view factor analysis of the 6 sides of the box. This was done analytically using the 
view factor formula for two right plates with a common edge. The method of Gebhart was used to 
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Table II Comparison of the temperature at the center of a sphere between the analytical 
solution and the numerical solution 


time 

(s) 

anal 

(5.0) 

anal 

(5.1) 

Ti- 

ll 

c 

n-6 

n= 11 

n= 13 








1 

99.6 

99.7 

93.3 

97.1 

99.4 

99.5 

2 

88.7 

90.0 

78.5 

84.7 

91.1 

91.2 

3 

68.9 

71.2 

62.8 

68.5 

75.1 

74.9 

4 

50.6 

53.2 

49.2 

53.2 

58.3 

57.9 

5 

36.5 

38.9 

38.0 

40.5 

44.2 

43.7 

6 

26.1 

28.2 

29.2 

30.6 

33.1 

32.6 

7 

18.6 

20.4 

22.4 

22.9 

24.6 

24.2 

8 

13.3 

14.7 

17.1 

17.2 

18.3 

17.9 

9 

9.5 

10.6 

13.1 

12.9 

13.6 

13.3 

10 

6.7 

7.7 

10.0 

9.6 

10.1 

9.8 

11 

4.8 

5.5 

mm 

mm 

mm 

mm 

12 

3.4 

4.0 

5.8 

mm 

5.6 

5.4 

13 

mm 

2.9 

4.4 

4.0 

4.1 

4.0 

14 


2.0 

BEM 

3.1 

3.1 

3.0 


mathematically compute the radiative interchange factors and the resultant K ra(1 matrix. 

Case 2: 6 Node, Discretized Radiative Transfer 

The internal radiative transfer was modeled by discretizing each face of the box. Ideally this 
will increase the accuracy of the radiative interchange. However, since the AutoSINDA radiative 
transfer algorithms only tabulate this data in the 6x6 matrix, this very large set of data had to be 
summed over the subelements within each face to back out the net radiative interchange for each ot 
the six faces. 


Radiative Transfer Validation Results 

AutoSINDA was executed on the 6 sided solid geometry model using 10 thousand rays per 
element, allowing a maximum of 30 reflections per ray. A mixed mode reflective model was chosen 
allowing both photon analog radiative transfer propagation as well as the continuous ray progression. 
The execution time for this case is on the order of 30 seconds wall clock time for a SunSPARC 10. 
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Table III Physical Model of the 3x4x5 cm Box for Validation 


Box 

Face 

Outward 

Normal 

Face Area, 
cm z 

Material 

Broadband 

Emissivity 

1 

-X 

20.0 

Cast Aluminum Alloy, Type 360 

0.200 

2 

+X 

20.0 

Aluminized Kapton 1, Kapton side out 

0.740 

3 

-Y 

15.0 

Carbon Black Tedlar Film 

0.800 

4 

+ Y 

15.0 

Highly Reflective Aluminum Alloy, 
2024-A1 

0.030 

5 

-Z 

12.0 

1 

White Paint, Chemglaze A276 

0.881 

6 

+ Z 

12.0 

Black Paint, Chemglaze A382 

0.975 


Table IV shows the results of AutoSINDA in comparison to both of the methods described 
above. Overall, the mean deviation for the six node view factor method is about 4% with a 
maximum deviation of 24%, so the comparison against the simple 6 node model is quite unacceptable. 

However, using the more accurate many node model, the mean deviation falls to only 2% with 
a maximum deviation on the order of 5%, an increase of almost 5 times the accuracy. Hence it is 
concluded that the AutoSINDA algorithms compare very well with the more accurate method of 
computing the radiative transfer. 


SAMPLE OF COMPLICATED GEOMETRY 

As a sample of a complicated geometry, Figure 2 shows a spiral antenna attached to a flat plate 
which has been subjected to orbital and laser loadings. This example is not presented to compare to 
analytical or measured data, but merely to show the range of applications and the usefulness of data 
visualization. 
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Table IV Results of the Radiative Transfer Validation 


Box 

Faces 

AutoSINDA Results 
[W/K 4 ] 

Analytic Results 

6 Surface Elements 

600 Surface Elements 

Krad 

%e 

Krad 

%t 

1-2 

7.06995xia 12 

7.22104xia 12 

2.14 

7. 16098x1 0 12 

1.29 

1-3 

5. 04044x1 O' 12 

5.10710xia 12 

1.32 

5.10193xl0 12 

1.22 

1-4 

1.65488x1 O’ 13 

1.68235xia 13 

1.66 

1.73681x10-" 

4.96 

1-5 

4.520l9xia 12 

4.55183xia 12 

0.70 

4.56767x1 O' 12 

1.05 

1-6 

4.87479X10- 12 


4.47 

5.06076xl0 12 

3.82 

2-3 

2.1I906xia" 

2.18800xia 11 

3.26 

2.17028x10" 

2.42 

2-4 

7.23058xia 13 

7.20755xia 13 

-0.32 

7.32650xl0 13 

1.33 

2-5 

1.901 17xia M 

1.9501 lxia 11 

2.58 

1.93365X10 11 

1.71 

2-6 

2. 09295x1 a 11 

2. 18181x1a 11 

4.25 

2.14510x1a 11 

2.49 

3-4 

5.39 125x1 a 13 

5.66541xia 13 

5.09 

5.40254xl0 13 

0.21 

3-5 

1.47630x1a 11 

1.56077x1a 11 

5.72 

1.53137x1a 11 

3.73 

3-6 

1.68516xia n 

1.74621xia u 

3.63 

1.69953xlO u 

0.85 

4-5 

5.11 129x la 13 

5.14137xia 13 

0.59 

5.19854xl0 13 

1.71 

4-6 

5.64379xia 13 

5.75225xia 13 

1.92 

5.75496x1 O' 13 

1.97 

5-6 

1.10285xia u 

1.36538xia n 

23.8 

1.12036x10" 

1.5 
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AUTOSINDA FLOW CHART 
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Figure 1. Flow Chart of the AutoSINDA Process. 
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